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ABSTRACT : 

In this paper we discuss three known methods for obtaining and 
updating the modified Cholesky decomposition (MCD) for the particular 
case of a covariance matrix when one is given only the original data. 
These methods are the standard method of forming the covariance matrix 
K then solving for the MCD, L & I) (where K = LDL 1 ); a method 
based on Householder reflections; and lastly, a method employing the 
composite-t algorithm developed by Fletcher and Powell ( Math Comp . , 

28, 1974, pp. 1067-1087). For many cases in the analysis of 
remotely sensed data, the composite-t method is the superior method 
despite the fact that it is the slowest one, since (1) the relative 
amount of time computing MCD's is often quite small, (2) the 
stability properties of it are the best of the three, and (3) it affords 
an efficient and numerically stable procedure for updating the MCD. 

The properties of these methods are discussed and FORTRAN programs 
implementing these algorithms are listed in an appendix. 

Institute for Computer .Services and Applications 
Rice University 
Houston, TX 77001 
May, 1976 
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Introduction 


In digital processing of remotely sensed data, as well as many 
other areas employing multivariate analysis, solutions to many of the 
problems are formulated in terms of covariance matrices. Often 
these solutions are expressed in terms of linear transformations 
involving a covariance matrix or its inverse. In these and other 
cases, it is often more sound computationally for one to employ the 
Cholesky or modified Choleskv decomposition of the covariance matrix 
rather than the original matrix itself^ 1 \ Even in cases where the 
original covariance matrix is to be modified by the addition or dele- 
tion of data, it still may be computationally prudent to utilize these 
decompositions. 

The purpose of this paper is to discuss methods for computing 
and updat* the modified Cholesky decomposition (MCI)) of a covari- 
ance matrix. These methods will be examined from the point of 
view of their ability to update the MCI) when data is to be added 
or deleted, as well as computational efficiency and numerical 
stability. In particular, three methods for accomplishing the above 
will be discussed: 

1) Standard--one computes the covariance matrix from 
the defining equations and then calculates the MCI) of it. 

2) Householder- -here one directly computes the MCI) 
of the covariance matrix from the data using Householder reflec- 
tions^. 

3) Composite-t--this method was developed by Fletcher 
and puwell^) from work previously done by Bennett^ and 
Gentleman^. The method essentially uses Givens rotations^) 0 f 
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the data one point at a time to directly compute the MCI) of the 
covariance matrix. Updating is straightforward and efficient. 

Table l summarizes the properties of each of these methods. 

Here n is the dimension of the data and m is the total number of 
data vectors. Though numerical stability may not play much of a 
role in most cases, the times when it d«>es, may occur without the 
user being aware of any difficulties. Thus this situation may lead to 
erroneous interpretations of the results. A method for computing the 
MCI) should be chosen with this in mi id. Also, in many areas of 
digital processing of remotely sensed data, the actual computation 
time for computing the MCI) is inconsequential, so that the composite-t 
method with its superior stability, may be optimal despite its relatively 
slow performance. The added benefit of an efficient and stable updating 
capability may al >o lie of value. 
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II. Methods for Computing and Updating the MCI) of the ("ova nance 
Matrix 

Given X the nxm data matrix containing m multi- 
variate n-dimensional vectors, the mean vector u is defined by 

m 

W j = V *ji J = 1 , 2 , . . . , n (I) 

1 = 1 

and the covariance matrix K, by 

K =m'rr I < x .,-u> (x.,-u)' (2) 

i = 1 

th 1* 

where x ^ . is the i l data vector and denotes transpose. 

K is symmetric and positive semi-definite. (It should also be 

noted that K is singular if men + I). The MCI) of K s 

T 

given by K = L DL where I) is diagonal with positive diagonal 
entries and L is unit lower triangular. 

A. Standard Method: 

The usual method for computing K comes from 

rewriting (2) using (1) to yield 

m mm 

k ij = m"-l [ I x i £ x j£ * rfi ( I x ijt)(I x j£ 

l = 1 £ = 1 £ = 1 

The MCD of this is then computed (sec c. g. ref. 6). This 

method requires (we consider cases where m > > n > 1) 

2 

approximately ■ ■ . - ■ multiplies to compute K anu another 

o ^ 

n /6 multiplies to compute the elements of L and D. 

Though K itself may be computed with acceptable precision, 
functions involving K amv be evaluated quite inaccurately since 



5 


matrix products of the form YY may be quite ill-conditioned' \ 
Updating L and I) in this manner is time consuming since one 

must first update K and then recompute L and D. Another 

method for updating L and D directly will be discussed in 
section C. 

B. Householder: 

One way to avoid problem of the possible ill 

conditioning of K is to compute L and I) directly from the 

data. This may be done by using Householder reflections on the 
data matrix as follows. 

Let M be the n x m matrix 



Then eq. (I) can be written 

K - jjtj (X-MXX-M)' 1 

If we then let 

X - M = R T Q (3) 

where R is m x n uoper triangular and Q is m x m and 
orthogonal, we may write 

1 T T IT 

K = R QQ R = R R 

m - I ' ' m - 1 

T 

- L DL 1 

where L and D are the MCI) of K as before with 
L - R 1 
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Rewriting (3), we have 


Q T (X-M) T 


iT 


R 


We may then write Q as a product of n Householder reflec- 
tions (see e. g re f . 2) 

,T 


Q 


P n P n-1 ••• P 1 


where P. annihilates all elements from i + 1 to 
n of the i tb column, changes the ii 1 * 1 element and does not 
change elements 1 to i-1. 

The algorithm, then for computing L and I) in this 
fashion is : 

Householder MCI) Algorithm 


m 


1. Compute u. = — £ 


i t 


t = 1 


i = 1 , 2 n 


2. Form t.j = Xj. - u { 


i = 1 , 2 n 

J = 1 , 2 , . . . , m 


3. For i=l,2,...,n, compute 

m y 

a) » = »gn (t,,) • (I t 2 ) 

j = i 


b) u " (0 ’ 0 C ii + a i» r i +1 , i, 

r m , i ^ 


c) e 


a u 
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m 

d) 'ii - i ii - ? u j I u k 'n 

k = t 

for J = i , i + 1 , . . . , m 

and for t - i + 1 , i +2 , . . . , n 

( N. B. denotes the replacement operation) 

e) t n - -or 

4. For i=l,2,...,n, compute 

a) d t = t^j / (m-1) 

b) For j=l,2,...,i-l, compute 

*u * ‘ji 1 'jj 

It should be noted that the elements of T, X, L, and D can occupy 
the same storage locations (though X will then be lost) and that steps 
3 and 4 can be combined. 

2 3 

This algorithm requires approximately mn - n_ 

3 

multiplications, which may be (for m much larger than n ) up to 
a factor of two times slower than the standard method. However, here 
the stability problems have been alleviated due to the use of orthogonal 
transformations. Storage considerations may be a problem with this 
algorithm, since it functions most efficiently only if the entire data 
matrix is in core. A sequential version of this algorithm may hi used 
(see Chapter 27 of ref. 7) which alleviates the storage requirements 
and provides for an updating capability, but at a cost of increased 
computation time. The next method (composite-t), however, yields a 
more efficient and stable algorithm. 
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C. Composite-t: 

This method is based on an algorithm developed by Fletcher 
and Powell* ** as a more numerical’ accurate extension of algorithms 
developed independently bv Bennett* 4 * and Gentleman*’’*. Essentially, 
Givens rotations*'** are used to directly compute the MCI) from the 
data as in the previous method. Instead of working with all of the 
original data at once as in the Householder method, this algorithm 
updates the MCI) as each data vector is processed. In this section, 
we will present two algorithms which employ the composite-t method 
to calculate the MCI) and update it. 

The generalized composite-t algorithm for a rank one update of 
L and I) of the form 

LDL 1 ♦— LDL 1 + -L zz 1 

l l 

for a positive semi-definite matrix where we assume t } t 0 , that 
the rank of the matrix never decreases, and that t ^ 0 only if 

D has full rank (i.e. Y has full rank), «s: 


Rank-o ne Composite-t Algorithm 


1 . If t j > 0 go to 5 


2. Solve Lv = z for v 

(i.e. v j = z j , v . = z { 


i = 2,3,...,n) 


i- 1 


' t v 

ij ) * 


j = 1 


3. For is! , 2 n, compute 

‘i+i ■ l i + v ? ' d i 
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4. If any t ^ + j > 0 , then 

a) set t j « et Jf where e Is the machine 
precision 


b) 

for 

i = n , n - 1 1 



l i * *1+1 * v f 1 d i 

For 

i* 

* 1 i 2 | i t | n 

a) 

v i 

= Zj 

b) 

If 

dj 4 0 go to substep c) 


0) 

if Vj ' 0 to to sub-substi p (4) 


(2) 

*1 + 1 = *i 


(3) 

go to substep k) 


(4) 

d j « / tj, = Z* / Vj 


(5) 

calculation complete 

c) 

If 

2 

t j >0 then t + j = t j + v j / 

d) 

*1 

= 'i+i 1 'i 

e) 

d l 

- dj », 

0 

If 

i = n then calculation is complete 

g) 

e i 

= (oj Vj / dj) / t J+1 

h) 

If 

a j >4 then 


1) 

Y i = *i 1 *i +1 


2) 

for j = i + I , i +2 , . . . , n 



XX = V j * j j + 0 j z. 
Zj - Zj - Vj £ j j 


3) 

l j j - XX 
go to substep k) 
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i) z * — z 0 * Vj f * j 

j) ^ ? j z ♦ 

k) Return to subs rep a) 


Note that only the last n - 1 components of i * . ( i { . * I ) and 
n - 1 + 1 components of z «, need be involved in these computations. 
Note also that the number of multiplications performed in this algo- 
rithm is data dependent. A detailed error analysis of this algorithm 
is given in ref. 3 showing this algorithm to be quite stable. 

To employ this algorithm, we must first rewrite eq. 2 
expressing the covariance matrix, K ( r f 1 ^ associated with the 
first r + 1 da'a vectors to that, associated with the first 

r data vectors : 

K<r + t> . K<r> + z <r> z <r) T 

where 



Given m 
K (r+!) 


data points, then the algorithm for computing the MCI) of 
is : 
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Composite-t MCI) Algorithm 

1. Set L*I n , D=0, and s*=x*| 

2. For r« 2 , 3 , . . . , m 

a) set t| « r 

b) for i=1,2,...,n, compute 

Zj « 8j / v' r - 1 - v r- 1 x Jr 

8 i “ 3 * * * * 8 i + x ir 
(N. B. For r=m Sj = m u ( ) 

c) use Rank-one Composite-t Algorithm to 
update L and I) (note that t j >0) 

3. For 1=1,2 n 

dj - dj / (m - 1 ) 

The number of multiplications involved in this algorithm when m > > n 
is approximately mn + 7m n. m square roots are also necessary. 

After L and D have been computed, data vectors may be 
added or deleted yielding a modified L and D by using the following 
algotithm 

Composite-t Update Algorithm 

1. Compute d* - d „ * ( m - 1 ) where rh 
is the net number of points used to compute 
L and D. If s is unavailable, it may 
be computed from s * = m u * 
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2. If n point is being added, set t j = m+ 1 , 

y * » and q * m+1 

3. If a point is being deleted, 

a) set t ] = - m , y = ~1 . 

and q = m - 1 
L) s ^ •— fi ^ * o i 

where a * is the data vector to be 
added ^r deleted 

4. Compute z * = s * / y-yo* 

If a point is to be added, 

set S ^ •- S * + iy + 

5. Use Rank -one Composite- 1 Algorithm 

6. Compute d * - d * / (q-1) and 

s et m = q 

2 

This algorithm requires approximately 3n /2 multiplications to 

o 

delete a data point and n to add a data point. 

III. Numerical Examples 

In this section, we present some numerical examples illustrating 
the properties of the algorithms discussed in the previous section. 
Listings of the programs used to implement these algorithms are 
given in the appendix. 

Using some 12 dimensional pseudo- random data of 100 points, 
we tested all three methods on an IBM 370/155. All algorithms 
yielded the same results to ~5 decimal places. The times for 
computing the MCD’s by each method are: standard method— 
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. 59 sec., Householder — .78 sec., and composite -t — 1.02 sec. Note 
that the order of these timings is as predicted, but due to differences 
in bookkeeping and other ope rations involved in each method, these 
timings do not follow the ratios of the number of multiplications in 
Table I. 

To test the stability of the three methods for computing the 
MCD, we used the data matrix 

( 1 - .999 -.001 

1 - .99 -.01 

1 - 1 . .001 



which generates an ill-conditioned covariance matrix. The resultant 
L's and D's for each method and the exact L and I) are given 
below (to six digits): (We use the subscripts E, S, HR, and CT 
for exact, standard, Householder, and composite-t methods, 
respectively. Also only the below diagonal elements of L are 


given, in the order i ^ , t ^ 

» £ 32 > 


L E = ( . 995505, 

1.00050, 

-. 185148) 

L s = (.995504, 

1 . 00050, 

-.180695) 

L, ir = (.995505, 

1.00050, 

-. 185112) 

L ct = (.995505, 

1.00050, 

-. 185147) 


D b = diag. (.666001, .405405 x 10 ’ 4 , . 444444 x 10 " 6 ) 

D s = diag. (.666001, .405539 x 10 " 4 , . 823690 x 10 ' 6 ) 

D }1R = diag. (.666000, .405427 x 10 ' 4 , . 444646 x 10 ' 6 ) 

D ct = diag. (.666000, .405405 x 10 ' 4 , . 444447 x 10 " 6 ) 
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To illustrate the effect of these rounding errors, we then solved the 
system 

Kb = L DL 1 b 

The computed b's (accurate to six digits), are given below 

b E = (3182965., -518130., -2665833. )* 

b s =(171 5970. , - 283491 . , - 1 433039. ) 1 

b I)R = (3181339., -517794., -2664543.) 1 

b (T = (3182936., -518124., -2665813.) 1 

Note that b^ is off by a factor of ~ 2 (mostly attributable to the 
computed value of d.^ ), whereas bj )R is accurate to 3 digits and 

b^.j. to 5 digits. 

We next tested the updating capability of the composite-t update 
algorithm. When data points are added, the algorithm yields the same 
results as if one started with all of the data points since, except for a 
few multiplications, the computations are equivalent. When data is 
deleted, however, the answers may differ, since the process of data 
deletion is intrinsically less stable^. The following example illus- 
trates this : 

Let 

X = / 1 - .999 -.001 0. 1 

1 - .99 -.01 0. 2 

1 -1. .001 .001 1 
(This is the same data as used above with the addition of the data 
vector ( 1 , 2 , 1 ) 1 . ) Using o = ( 1 , 2 , 1 ) 1 and specifying deletion 
to the composite-t update algorithm yielded 
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L = (.995504, 1.00050, -.186673) 

D = diag. ( . 665^99, . 402583 x 10 " 4 , . 431288 x 10 " 6 ) 

which is to be compared to L ( .j and -| as computed above: 

L ct = (.995505, 1.00050, -.185.47) 

D ( , t = diag. (.666000, . 405405 x 10 " 4 , . 444447 x 10 ' 6 ) 

The differences are in the second and third digits of some of the 
computed quantities. It should be pointed out, however, that this is a 
particularly ill-conditioned example, and other examples yielded satis- 
factory results. 

IV. Conclusions 

Though efficient, the standard method suffers from an 
inability to update accurately and efficiently the MCI) , as 
well as stability proolems associated with having to work with 

T 

matrices of the form YY . The Householder method obviates 
these problems at the cost of storage requirements and efficiency. 
Though slower still, the composite-t method drastically reduces 
the storage requirements, readily provides for updating of the 
MCI) and improves computational performance from a stability 
standpoint. Which method one should use depends on the problem 
at hand and the weights one assigns to the various trade-offs 
between speed, stability, and updating capability. 

In many of the computations for the analysis of remotely 
sensed data, the actual calculation of the covariance matrices 
and their MCD's takes relatively little time, so speed may not 
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be an important factor. In this case, the optimal choice would appear 
to be the compostte-t met^'Kl, due to its superior numerical stability, 
relatively small storage requirements, and its updarinc capability. 

In areas such as signature extension, the updating capability of this 
method could be especially valuable. 
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APPENDIX 

Listings of the program used to test the three methods are given 
below. KBYSM computes the MCD by the standard method (subroutine 
MCHLSK is used to actually compute the MCD from the computed 
covariance matrix). KBYIIR computes the MCD by the Householder 
method. KBYCT computes the MCD by the Composite-t method, 
usirg subroutine COMPT which computes a rank one update of the 
MCD (Note that all of the data is in KBYCT, though the algorithm 
only requires that one data vector at a time be available. 'Phis was 
done for timing purposes only.) CTUPDT updates the MCD using 
subroutine COMPT. 



r.r>r o r,m non nrnrftftOftftfiftn nnnn 


18 


SLERLUT INE KHYSM(Ai4 t N«MAN«LDI 

THIS ROUTINE COMPUTES 1HL MCD UE A LUVAi' | A NCI MATRIX HY T HT 
S T AND AMO ML THUD 

ME AL • A X ( M XN * 1 I * L D ( I ) 

RF ALArt SI ( I 2 ) . S2 ( 78 ) 

A - THt N DV M DATA MATRIX WHOSE F |M:,T DIMENSION IN T-1 1 C ALL I NO 
PROGRAM |b MX N 

M • I Ft NUMUEW OF DAI A vectors 
N - THE DIMENSION OF THE DATA 

LL) - f HL RESULTING MCU CONTAlNGG THE FLLMLNTG OF L 6 l)T THIS 
MATRIX IS ST UHL O IN SYMMETRIC STORAGL MUOt (l.t. I Q*ER 
TRIANGULAR PORTION STORED bV ROWS) * I Til THL ELFMf NTS OF C 
OCCUPYING THF DIAGCNAL ENTRIES* 

INITIALIZE 

DO 20 1*1.!* 

20 bill 1-0. DO 
UU JO I - I « 7 tt 
JO S2 ( l l>U.OC 

COMPUTE THL SUM UF T Hi. l.ATA VECTORS AND I Hr I R C *>OSS-PROOOC T S . 

DO 1C L-l.M 
K - C 

DO : 0 1= 1 *N 

SI (1 )*S1 ( I ) AX( I *L ) 

00 10 J=1.I 
K*K* I 

10 S2(K)*S2(K)+X( I i L) ♦*( J (L I 

COMPUTE THL COVARIANCE MATRIX G STURL IT IN LO* 

K= C 

DO AO 1=1 .N 
DO 4 0 J* 1 • I 
N * K ♦ 1 

1.0 (N) =( S2< EO-S1 ( 1 ) *bl ( J )/M )/( M- I ) 

40 CUMI NUt 

MCHLSK COMPUTES THE MCL UT THE COVARIANCE MATRIX G OVrRW^ITLS 
THE RESULT UN IT. 

CALL MCHLSMLD.N.S1 .S2 ) 

RE TURN 
END 


. iOIMAL page is 

OF POOH QUALITY 
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SUEMLUllNt »*CHLin(»U,NV.i)U^.Ot I) 


iMIi UUUIlNfc CiiMPull j TK MUDIFILI) LMULtS>V UC CLMPO » I T I ON l.» 
f HL CUVAHI anil matrix. IHL DFCLMPObt T It Kb . VfPLAr !H» fLt^'MS 
ur r hi cuva.< i ANct matrix, 

KK.iL O L • 


kk - rnt cuvawianci •'Air lx stoned ik syMM^r^ic su dagf vet*. 

NV - T ME NUMufcH U» LhANKILa Lt-fO 

OUM - A UOUULt PHECIGIOK aOMK AH* A <JF A I 7F NV-I 
Otf - f Hi ULUMMINAM O F Thf CUVAWIANC^ K A T P I X » 

k£ AL KKI | ) 

It Al Art OUM ( | ) 

ME AL H.k I . U . Tf 
LOGIC AL • 1 JLI 
JL 1= . I Mil . 

J l *0 
JO*0 
OE 1=1# 

LOUP OVEH ALL CHANNlLb 

l)U 10 J 3 l » K V 
KL = J- I 
L = J ♦ J 
JO=Jl 
J 1 = J 1 ♦ J 

It =0 .oc 

I F <Jtl) GO TO 12 
K l =0 

CCMFUTF ? Hh DIAGONAL tL»MtMS OF l> AND r jlli*t IN Kc 
ItMPimAti ILY slUNl. IMr- UfUDUCI KK ( I . I I ♦ K < < J , | » IN UuMII 


DO lb 1 = 1. KL 
M = KK < JL ♦ I I 
K I =K I ♦ I 
M | =KK ( K I I *K 
rt *TF -W 1 *P 
DUM( I ) = M| 

15 CUM INUL 
12 CUM INUE 

1F = TFFKK< J| ) 

K K ( J | » = TF 

DE 1= Ufc 1*TT 

IF (L.GT.NV) gu TU 10 
I M C = J 1 “L ♦ I 

COMPUTE T HL H.J-TM CLEMENT Or L UL> I NO I I 


DU 20 IH-L.KV 
I MC= I PU* I H- | 

T I =0.00 

IF (JLI) GU Hi I f. 

DU 2b I =1 . KL 
T l = T 1 -DUM< I ) ♦ K K ( I H O ♦ I ) 

2b COKi iNLt 

1(3 KK ( I HU* J ) = ( T l FKK ( 1 GD4 J ) )/T f 
20 CUNT INUF 

Jfc 1= .F ALSt . 

10 LUMlNUc 
HE TUKN 
eND 
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St- EPUG i INI AUYHP(X«WO«N*MtL.MXN*LOI 
Ht Ac ♦ A LO( I ) 

PlAC»A A ( MAN ( M ) , MU ( N ) • U ( 4 I 
PtACAH S 
C 

c ui vi n ihl matpix of uhslpv at il ns a. cumpuii tmi m f an a <o mc n c r tmi 

C Pt SOLTI No CJVAP|ANC( ma I ■< I * f. sluPI IT IN f HI TPIAnC.UIAP 

C HAST OF A. HUUbtHULJLh <LM.rtUUN'j Axl USl • T( CCMPUlc Th ( M(IJ. 

c 

C X • DATA MATPIX WHICH lb Ul blt-LYf C 

C MU - CK 00 T HU i • f H l M| AN VH.I0W 

C im - I He OIWtNSION 

C M • THC NUML'tP Uk UMbtPVATICNS TU ML U SI L> 

c MAN - I ML OlMtNSlGN ( • ok MUWj ) (If X IN I HI CALLING PPIIG. 

C U - wUNkIng STUpAGL OF OIMLNSICN AT ((AST v 

C CO - O t S I ml L MCO sTOntD IN SVM STUPAGL Moot 

c 

NP i = M- | 

C 

C SStT CP A A TP IX T C Ok TP I ANOUL A«|/tn 

c 

Ou 10 I * 1 • . 

s*c.uo 

c 

C COMPOTE M.ANS 

c 

00 4*0 J = 1 . M 
20 S= S» A < I . J I 

MO < I > = £/(• 

c 

C COMPOT T Ht MATHlA X-MkANS 

c 

OU JO J=I.M 
JO X( I • J I = X ( 1 i J l-MU( | ) 

10 CONTI NOt 

c 

C PEPkOPM HOC SLMULOl h T M A NSk LPM A T I LNS 

xx s 0 
C 

OU AO I 3 1 • N 

c 

C COMPUft NtCkSSAwy oUANTITIkS TO ANMHlLATc l»u T T oM PACT Ok l-TM COL 

C 

S» C . O 0 

c 

C ONLY LAST M-Ik| tLLMtNTu ilk U A PL USLO 

C 

OU AG J = I • M 
AX = A ( I • J I 
£♦ X X AXX 
A 5 O ( J I = A A 

ALP-^IOM jNOLI OS OPT ( s I I .U< III 
It ( I . E C • W I GO TO A A 
U( 11-0(11 k ALP 
ML TA= AlP*UI I I 

1 I «l ♦ 1 
C 

C APPLY TWANSFUPMAT ICN TO HOWS I ♦ A To M (, COLS Ik I TO N 6 SIT 

C l.l-TH tLkMtNT TO -ALP 

C 

DC bC L* I I • N 
S*0.D0 
Oo bb K= I . M 
65 s= Sk'U ( x I #X( L »x) 

XX=b/BfcT A 
OU 60 J— l • M 
X(LiJl : X(LiJ)-U( J I ♦ A A 
50 CL NT 1 NOt 
A A A « I • I 1 — • AL P 

Ik ( I .tU. I I GU TU A<? 

c 

C COMPUTE LCD 'POM L* SOP I ( 0 ) * W-TP ANS/ S UP T < M» l I . STUPT L IN 

C LUWWP TPIANkULAP PAPT Lk A t O AlUNo O I A(. . 

C 

11=1-1 

OU t.C J= l . 1 l 
KK =XX k 1 

bC LO (xX )■= x ( I • J l/X ( J . J I 

4<! KK. -XX ♦ 1 

LU(XX)=ALP*ALP/NP| 

AO CU NT 1 NUt 
PE TUPN 
k n n 
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SutHuUTINl CUMFT dOd.^tNi V.TFP) 

Wt*LM Lid ) 

kl*LM I I I I • l I N I . V ( N ) 

KtOL*8 fFPIN) 

HEAL ♦ S 

LUC IL AL * I ILLS > WNUtHH tL ALP 

c 

C TH 1 1* HUUI INI IS AN IMPl FNtNTAT ION OF T Ht CUMP.SlTt - I ALGCMTHH 

C TU F*H RF uH F A HANK I UPDAIT OF T FI WCU sTUWtl* IN AHFAY LDI I.C. 

C K*L*L*L-THANi t «ft *I3M TU COPPUffc L • LO • S • T .K • «K4 i * 7-TH AK £ /T I | » 

C C K**L**U* *L * •! k ANSI • 

C LD - ARRAY CONTAINING L C C SlUHFO IN SVM.SlUHAol MODI 

C I • AN N ♦ l VtCTUH KMUSt f' I H S f e Lk Ht NT IS AS AU;VF 

C l - VcCTuW LF THt LPDAt t AS At‘LVF 

C N • I Ft O I FLNbluN 

C V - WiJHKlNG STU.-AUt UF LtNOTH . l,t. N 

C IMF - OULMLt PHLCISlON aORKINU SIOHAUF OF LINGTU .C.F. N 

c 

t pcs= r < i ) .or .o* 
if ( tpcsi cu in Jb 
c 

C A PUIM IS TU Uf Dtu no 

c 

EPS*S.Wfc-fe 

c 

C SUL Vt L*V±2 F UH V 

c 

KS 1 

V( II = /< 1 I 
UU 10 I»<>|N 

I J* I- I 

S* C.UO 

UU lb JB 1 . I J 
k = r.F I 

15 S C S4LU(K)*V(J) 

K = K ♦ 1 

VI I ) = Z ( I )-S 
10 CUNTINUt 

c 

C CU FPU I E IFfc T<|*S) 

C 

Kb C 

RNUfcl-H=.FAL St . 

DU £0 I * 1 « N 

KbKF I 

fMFl I )BV( | ) * V ( I I /L0<N > 

T < I F 1 ) = T I I )F1M»M I ) 

IF ( T ( I ♦ 1 ) ,tt .0. > KNDtRK-. TMJI . 

HO CUNTINUt 

IF I •KCT«KKDLKHJ UU TU .1 > 

C 

C WUUNDINU F.FNUW I • A S MAl/t A T(|Fl).Ct.O. SO C 1>«L C T I OH THIS 

C 

T I NF I ) =EPS4 T C l ) 

UU JO J = l * N 
I sN-JFI 

T( I) = T( IFl l-TMMl I ) 

JO CUNT I NLt 
35 COM iNLt 
I J B 0 

DU AC ! b 1 | N 

ll-UI 
I J= I JF 1 
V< I I - / ( I ) 

U I = LU( 1 J ) 

IF (Dl.ol.O. I Ou TO A A 
C 

C Dll) BO. su FANK CF D kill L I T F fc W I NC FT A St ON 1 I MAIN UNCO ANCF n . 

c 

IF ( UABSI V C 1 ) ) .01 * I .F-JC ) GU TL Id 
C 

C HANK OF C KILL H f MAIN UNCHANULL 

C 

T ( IF I ) bT< I ) 

GO TU AO 



rn nr n nor nnr 


22 


MAKK CF 0 » ILL INCMAUL t»V I 

42 LUIIJI>V(I >*V< I I '» ( I) 

IF (I «LU>N ) HE rUMN 
A* 1J 

UU 4b J: I I iN 
KiRAjal 

LUIkM^I J >/V< I I 
4b CuMIKCfc 
HE I OH N 
44 CU M 1 NUt 

UP C AT t 0 

IF llPCbl I I 1 ♦ I > * I (I I ♦ V ( I > *Vl I I/O I 
ALP* 1 ( | 41 1/ l( I I 
LOIUI -U I ♦ Al I’ 

IF ( I .K.M Ht fUMN 

UP LA T E L C MOU1FV l 4CCLH0IM.LV 

Ut T A = < V « I )/L I )/T ( I F| » 

LALP*.FALbt . 

IF ( AUM.LE.4i I OU 1C bif 

Thli» Ht 1 FUL CotO 1C INSUHL jTAIIlltV IF AL PPA 01. 4 

(.ALPS • IHUt • 

CA P> 1(11/1(1 + 1 I 
K* 1 J 

OU 1)0 J* I I »N 
K * A ♦ J* I 

XXIOAMILUI A) FUI IA«/I j| 

/ ( J ) -L ( J )- V ( ll»LU(KI 
LU(K)>XX 
t> C CUM l ACt 
GJ 1U 40 
a* I J 

DU t>0 J* 1 I » A 
K-* A * J • 1 

t ( j) = / ( J I- V ( I ) *U ( K ) 

LD ( A ) =LL ( A) FUt !.*•/( J) 

(lO CuM INUF 
4 0 COMIM? 

Hfc 1UI< N 
FAC 
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SOOHOOTlNt AUVtl I LO.N. van .< . X . N . T I 
Ht AL • 4 Ll)( I liblNIiK^kNiMl 
nt Ac** Mil 
INlUtlkU F 

IH1* HUJlIM COMl'oltj Iml MCI) OF THE ( l,VAt< I ANC*. MATRIX (, Mf AM 
> u I V t N I Ml N-Ol MtNSlwNAL M DATA VK1u»S jtlit-MI IN X. A_ C I ) AM 

STCWLO IN LL & IMF ANj IN VLCTllH S 

LO - LN UUTPoT. T Ht MCD OF T Ht 10VAH, MAMw|X STlMF'O IN 3V M . 

3 T OHAul FoOL » I TH 0 A«_ONO IH 0 I AoONAL • 

N • I Ht UMlNiluN 

MaN • flit f.uVMLW OF HCJA j OF X AS D|M| NSllM 1/ IN T HI l AH I No PF Uf. • 
S - IN uoTPOt. THt VlUM CF MF ANS 

X - I Ht DA I A MATH I X 

M - THE NUMutM OF OATA U.ttOFS TO Ot OSFO. 

T • M O H K I N Ci S TO H AO F OF DIMENSION .OL • A*N»l 


IMUAll/l LOU MATHltLS C VtClOH S 

IJ»0 

00 10 l*I.N 
I J* I J ♦ I 

LO ( I J )«0 

S( 1 ) = X ( 1 • I I 

IF ( I «t0. I I OO TO 10 

1 l •* 1“ I 

OO IS 0*1.11 
IS LO ( I J- I * J 1*0 • 

1 0 CONI I NOL 

LOOM U VI N ALL POINTS TO CJMFoTI L t O F OF (F-IJAK 

00 dQ R*2.F 

1 ( I > = FJ 

SK 1 * 3 0 F I ( F L O A I ( H — 1 )i 

COPPOTt l F On THIS x I. OPOATt t 
l)u dt> 1 = 1 • N 

T ( NF !♦ I )-3( I )/SH l-swi *X ( I iM 
d b S ( I ) = S ( I 1 ♦ X ( I * H ) 

UP CATE L L D 

CALL LCPPT (LO. T • T ( N»^ > .N. T ( i' *K ♦ I . T ( j • ? > ) 

dO COM INOt 

MUCH Y 0 S.T. K=L ♦OFL-l HAN £ t SI OH I MEAN IN S 

SR |*M» 1 
I J =0 

00 JC I * 1 • N 

1 J*1 JF I 

OO I I J I -LU ( l J 1 /SW I 
S ( I I = S ( I )/N 
JO CoNTINOt: 

Ht TUHN 
t NO 
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C I UPDT (lO.N.'A, ALL. b.ALP.T I 
NtAL*4 Lit I l.AHMM. 'j|M 

Ml ADH HI) 

LUCK AL * I Aull 

Tnlh MLIIIM UPL«lt> THL MCL bTLKCO IN LU MV I I T » ♦»■ - AU'TINl 
UM uLLtllM. A L)m I A VFCTOU Ah 'jl'l C It tLi> *> T H| LtGlCAl VALIA L( A. >C 

LU • IhL MCL bTOHLL IN SVM.jTOMAGt Mljfif WITH I) Ai'AL t Hi DlAC.tNAl 
N - Ttt U I Pt Kj I UN ot T Mi. LATA VtCTLM 

M - let Ml NUMulM Ut 0 A I A Ct.CH.Mh uSI U lu CuMCUlt CO* N t ILI 11 
UP L A ILL IN OUTPUT 

ADC - *1 It LAI » H. Mf Al U/ I L . «f |t DATA |( Cl Util 110. 

& • VtClLL UNTAININu t • ML ANL* 

b * ILL U Ut'UATLO UPCN tflOLN 
ALP - LATA VL C TUN IU ttf AO LI U/ UL LL T I D 
T • »UHK I NO UVUUAbt lit UlMiNhllN »GL.4*N* I 

MUOIt V U b • T • |M>| )*ML*OtL«THAKb 

XM sM 
I J *0 

00 10 1*1 »t 
I J- I J* I 

LU I I J I *LU ( I J I •< M- | ) 

10 LuMISUt 

It ( .NC t .ALL I l*L TO 1/ 

AOC A FCINT 

T < l UM ♦ 1 
V* SUM T < XM I 
0*Mt 1 
Uli TU 14 

DM. Lit A PLINT 

12 T I I ) --► 

V = SUL T ( XM- 1 » 

U-t-l 

CUMPUT1 / AND UPUATr b 
14 DU 20 1= I .K 

It (.Ml. AIL) S( I )»S( I ) -ALtl I I 
T ( ItNt I ) = b( I l/V-V • At_P( I ) 

It ( ADD b( I >-b( 1 ) «ALP( I ) 

20 CUNT I N0'_ 

UPDATE L L L 

CALL CCMt-T (LU. T • T ( Ntc ) .N. T I ? • N ♦ #* ) » T ( 3 *N ♦ i II 

MOU II V C b.T ,N*L • L • L — l K AN j 

1 J»0 

(J J 30 I- 1 . N 

I J* I J ♦ I 

LU ( I J )*LU( 1 J )/( l- I . ) 

30 CONTINUE 

Ul»I.T M 

Ms C 

ML lUKN 
t NC 
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